library(SpatialEpi)
## Loading required package: sp
library(readxl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
IMR1=read_excel("~/2018 Reported Infant Death.xlsx",skip = 1)
Malesex=read_excel("~/2018 Reported Infant Death.xlsx",sheet = "Gender Male")
head(IMR1)
## # A tibble: 6 × 5
## Districts `Place of recorded death` Cases Population (Total Infant B…¹ Age
## <chr> <chr> <dbl> <dbl> <chr>
## 1 Gaborone Gaborone 57 5986 >28 …
## 2 Gaborone Gaborone 23 2416 1mon…
## 3 Francistown Francistown 28 3678 >28 …
## 4 Francistown Francistown 19 2495 1mon…
## 5 Lobatse Lobatse 8 496 >28 …
## 6 Lobatse Lobatse 10 621 1mon…
## # … with abbreviated variable name ¹​`Population (Total Infant Born)`
d=aggregate(x=IMR1$Cases,by=list(Districts=IMR1$Districts),FUN=sum)
M=aggregate(x=Malesex$`Gender (M)`,by=list(Districts=Malesex$Districts),
FUN=sum)
head(d)
## Districts x
## 1 Central 303
## 2 Chobe 6
## 3 Francistown 47
## 4 Gaborone 80
## 5 Ghanzi 22
## 6 Jwaneng 6
names(d)=c("id","Y")
names(M)=c("Districts","Male_sex")
pop=IMR1$`Population (Total Infant Born)`
cases=IMR1$Cases
n.strata=3.25
E=expected(pop,cases,n.strata)
d$E=E[match(d$id,unique(IMR1$Districts))]
d=merge(d,M,by.x="id",by.y="Districts")
d$SMR=d$Y/d$E
head(d)
## id Y E Male_sex SMR
## 1 Central 303 24.67221 158 12.2810226
## 2 Chobe 6 48.55638 3 0.1235677
## 3 Francistown 47 52.64422 25 0.8927855
## 4 Gaborone 80 184.15797 42 0.4344097
## 5 Ghanzi 22 18.08976 12 1.2161576
## 6 Jwaneng 6 26.18266 4 0.2291593
Districts = raster::getData("GADM", country = "Botswana", level = 1)
## Warning in raster::getData("GADM", country = "Botswana", level = 1): getData will be removed in a future version of raster
## . Please use the geodata package instead
library(sp)
rown=sapply(slot(Districts, "polygons"), function(x) slot(x, "ID"))
rownames(d)=rown
identical(rownames(d), getSpPPolygonsIDSlots(Districts))
## Warning: use *apply and slot directly
## [1] TRUE
map= SpatialPolygonsDataFrame(Districts,d,match.ID = TRUE)
head(map@data)
## id Y E Male_sex SMR
## 1 Central 303 24.67221 158 12.2810226
## 9 Chobe 6 48.55638 3 0.1235677
## 10 Francistown 47 52.64422 25 0.8927855
## 11 Gaborone 80 184.15797 42 0.4344097
## 12 Ghanzi 22 18.08976 12 1.2161576
## 13 Jwaneng 6 26.18266 4 0.2291593
library(leaflet)
l=leaflet(map)%>%
addTiles()
pal=colorNumeric(palette="YlOrRd",domain=map$SMR)
l %>% addPolygons(color="grey",weight=1,fillColor=~pal(SMR),
fillOpacity=0.5) %>%
addLegend(pal=pal,values=~SMR,opacity=0.5,title="SMR",
position="bottomright")
labels=sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/>
Male sex proportion: %s <br/> SMR: %s",
map$id,map$Y,round(map$E,2),map$Male_sex,round(map$SMR,2)) %>%
lapply(htmltools::HTML)
l %>% addPolygons(color="grey",weight=1,fillColor=~pal(SMR),fillOpacity=0.5,
highlightOptions=highlightOptions(weight=4),label=labels,
labelOptions=labelOptions(style=list("font-weight"="normal",
padding="3px 8px"),
textsize="15px",
direction="auto")) %>%
addLegend(pal=pal,values=~SMR,opacity=0.5,title="SMR",
position="bottomright")
library(spdep)
## Loading required package: spData
## Loading required package: sf
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(INLA)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loading required package: foreach
##
## Attaching package: 'foreach'
##
## The following objects are masked from 'package:purrr':
##
## accumulate, when
##
## Loading required package: parallel
## This is INLA_22.12.16 built 2022-12-23 13:24:10 UTC.
## - See www.r-inla.org/contact-us for how to get help.
nb=poly2nb(map)
head(nb)
## [[1]]
## [1] 2 5 8 9 11 12 13 16
##
## [[2]]
## [1] 1 12
##
## [[3]]
## [1] 11
##
## [[4]]
## [1] 8 9 14
##
## [[5]]
## [1] 1 7 9 12
##
## [[6]]
## [1] 15
nb2INLA("map.adj",nb)
g=inla.read.graph("map.adj")
map$re_u=1:nrow(map@data)
map$re_v=1:nrow(map@data)
formula=Y~Male_sex+f(re_u,model="besag",graph=g)+f(re_v,model="iid")
res=inla(formula,family="poisson",data=map@data,E=E,
control.predictor=list(compute=T))
summary(res)
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 0.854, Running = 0.255, Post = 0.15, Total = 1.26
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.776 0.207 -1.190 -0.775 -0.369 -0.772 0
## Male_sex 0.020 0.004 0.012 0.020 0.029 0.020 0
##
## Random effects:
## Name Model
## re_u Besags ICAR model
## re_v IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for re_u 18852.32 18808.97 1290.87 13180.04 69010.69 3551.04
## Precision for re_v 3.03 1.26 1.21 2.82 6.10 2.42
##
## Marginal log-Likelihood: -93.87
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
head(res$summary.fitted.values)
## mean sd 0.025quant 0.5quant 0.975quant
## fitted.Predictor.01 12.3202871 0.70826119 10.98960406 12.2999794 13.7666008
## fitted.Predictor.02 0.1859291 0.05803335 0.09468564 0.1784237 0.3202189
## fitted.Predictor.03 0.8854982 0.12658889 0.66326213 0.8765663 1.1586771
## fitted.Predictor.04 0.4495950 0.04902685 0.36097573 0.4469962 0.5530276
## fitted.Predictor.05 1.1133586 0.23728791 0.72087100 1.0886917 1.6468377
## fitted.Predictor.06 0.2946768 0.09531422 0.14812262 0.2813314 0.5180600
## mode
## fitted.Predictor.01 12.2594222
## fitted.Predictor.02 0.1640752
## fitted.Predictor.03 0.8589968
## fitted.Predictor.04 0.4418447
## fitted.Predictor.05 1.0410090
## fitted.Predictor.06 0.2564490
map$RR=res$summary.fitted.values[,"mean"]
map$LL=res$summary.fitted.values[,"0.025quant"]
map$UL=res$summary.fitted.values[,"0.975quant"]
pal=colorNumeric(palette = "YlOrRd",domain=map$RR)
labels=sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/>
Male sex proportion: %s <br/> SMR: %s <br/> RR: %s (%s, %s)",
map$id, map$Y, round(map$E, 2), map$Male_sex, round(map$SMR, 2),
round(map$RR, 2), round(map$LL, 2), round(map$UL, 2)) %>%
lapply(htmltools::HTML)
leaflet(map) %>% addTiles() %>%
addPolygons(color = "grey", weight = 1, fillColor = ~pal(RR), fillOpacity = 0.5,
highlightOptions = highlightOptions(weight = 4), label = labels,
labelOptions = labelOptions(style = list("font-weight" = "normal",
padding = "3px 8px"),
textsize = "15px", direction = "auto")) %>%
addLegend(pal = pal, values = ~RR, opacity = 0.5, title = "RR",
position = "bottomright")
pal=colorNumeric(palette = "YlOrRd", domain = map$SMR)
leaflet(map) %>% addTiles() %>%
addPolygons(color = "grey", weight = 1, fillColor = ~pal(SMR), fillOpacity = 0.5,
highlightOptions = highlightOptions(weight = 4), label = labels,
labelOptions = labelOptions(style = list("font-weight" = "normal",
padding = "3px 8px"),
textsize = "15px", direction = "auto")) %>%
addLegend(pal = pal, values = ~SMR, opacity = 0.5, title = "SMR",
position = "bottomright")